"""
-*- coding: utf-8 -*-
@Time    : 2022/9/25 15:18
@Author  : issca_new
@File    : gnss_tools.py
@Software: PyCharm
@Attention:美色迷人心智，伤人根骨,瓦解人的意志
"""
import math
from math import *
from Tools import base_tools as bt
from numpy import *
from matplotlib import pyplot as plt
from Tools import mathtools as mt
import matplotlib.ticker as ticker
import matplotlib.colors as colors
import time


def isin(value1, value2, Ndata) -> bool:
    """
    由于卫星广播星历更新时间为2个小时所以n文件内大部分数据冗余，此处判断若同一历元，同一卫星的星历数据已经存入，则不再存入
    """
    for i in range(len(Ndata)):
        if Ndata[i]['卫星PRN号'] == value1:
            if Ndata[i]['历元'] == value2:
                return True
    return False


def head_num(lines_n:list) -> int:
    """
    获取头文件结束标志文件
    """
    # target_num :int= 0
    # for i in range(int(len(lines_n) / 3)):
    #     if lines_n[i].find('END OF HEADER') != -1:
    #         target_num = i
    #     else:
    #         break
    # return target_num + 1
    target_num :int= 0
    for i in range(int(len(lines_n) / 3)):
        if "END OF HEADER" == lines_n[i][60:73]:
            target_num=i+1
    return target_num



@bt.timeshow
def donfile(path) -> dict:
    """
    读取n文件
    """
    lines_n = open(path, 'r').readlines()
    head_nums = head_num(lines_n)  # 头文件结束标准
    Ndata = []
    data_line_num = int((len(lines_n) - head_nums) / 8)  # 总行数减去头文件/8，即历元数
    # 导航文件数据类型标签
    Ndataitemkeys = ['IODE', 'Crs', 'Delta_n', 'M0', 'Cuc', 'e', 'Cus', 'sqrt_A', 'TEO', 'Cic', 'OMEGA_A0', 'Cis', 'i0',
                     'Crc', 'omega', 'OMEGA_DOT', 'IDOT', 'L2Codes', 'GPSWEEK', 'L2PCODE', '卫星精度', '卫星健康状态', 'TGD',
                     'IODC', '电文发送时刻', '拟合区间', '备用1', '备用2']
    # 开始读取数据
    for i in range(data_line_num):
        Ndataitem = {}  # 字典，准备接收某卫星数据
        k = 0
        # 按8行读取
        for j in range(0, 8):
            data_countent = lines_n[head_nums + 8 * i + j]
            if data_countent[0] != "R":
                Ndataitem['数据组号'] = i + 1
                if j == 0:
                    Ndataitem['卫星PRN号'] = data_countent.split('\n')[0:3][0][0:3].replace(" ", '')
                    Ndataitem['历元'] = data_countent.strip("\n")[3:22].replace("  ", ' ')
                    # a = gpstime(julday(bt.all_float_list(bt.separate(Ndataitem['历元']))))

                    if isin(Ndataitem['卫星PRN号'], Ndataitem['历元'], Ndata):
                        # 由于卫星广播星历更新时间为2个小时所以n文件内大部分数据冗余，此处判断若同一历元，同一卫星的星历数据已经存入，则不再存入
                        break
                    Ndataitem['卫星钟差（s）'] = float(
                        str(data_countent.strip('\n')[23:42]).replace(" ", '').replace("D", "e"))
                    Ndataitem['卫星钟漂（s/s）'] = float(
                        str(data_countent.strip("\n")[42:41 + 20]).replace(" ", '').replace("D", "e"))
                    Ndataitem['卫星钟漂速度参数（s/s/s）'] = float(
                        str(data_countent.strip("\n")[61:]).replace(" ", '').replace("D", "e"))
                else:
                    try:
                        Ndataitem[Ndataitemkeys[k]] = float(
                            str(data_countent.strip('\n'))[3:23].replace(" ", '').replace("D", "e"))
                    except:
                        Ndataitem[Ndataitemkeys[k]] = 0
                    try:
                        Ndataitem[Ndataitemkeys[k + 1]] = float(
                            str(data_countent.strip('\n'))[23:42].replace(" ", '').replace("D", "e"))
                    except:
                        Ndataitem[Ndataitemkeys[k + 1]] = 0
                    try:
                        Ndataitem[Ndataitemkeys[k + 2]] = float(
                            str(data_countent.strip('\n'))[42:61].replace(" ", '').replace("D", "e"))
                    except:
                        Ndataitem[Ndataitemkeys[k + 2]] = 0
                    try:
                        Ndataitem[Ndataitemkeys[k + 3]] = float(
                            str(data_countent.strip('\n'))[61:].replace(" ", '').replace("D", "e"))
                    except:
                        Ndataitem[Ndataitemkeys[k + 3]] = 0
                    k = k + 4
            if len(Ndataitem) > 4:
                Ndata.append(Ndataitem)
    # 返回字典
    return Ndata


@bt.timeshow
def readnfile(path):
    f = open(path, "r").readlines()
    head_nums = head_num(f)
    content = f[head_nums:]
    del f, head_nums, path
    Ndataitemkeys = ['IODE', 'Crs', 'Delta_n', 'M0', 'Cuc', 'e', 'Cus', 'sqrt_A', 'TEO', 'Cic', 'OMEGA_A0', 'Cis', 'i0',
                     'Crc', 'omega', 'OMEGA_DOT', 'IDOT', 'L2Codes', 'GPSWEEK', 'L2PCODE', '卫星精度', '卫星健康状态', 'TGD',
                     'IODC', '电文发送时刻', '拟合区间', '备用1', '备用2', '备用1', '备用2', '备用1', '备用2', ]
    glonassitemkeys = ['x', 'xv', 'xa', 'health', 'y', 'yv', 'ya',
                       'frenum', 'z', 'zv', 'za', 'ageoper', '备用2', '备用1', '备用2', '备用1', '备用2']
    res, mid = [], []
    for i in range(0, len(content)):
        if content[i][0] in ["G", "E", "R", "C", "J", "I", "S"]:
            res.append(mid)
            mid = []
        mid.append(content[i][:-1])
    res.append(mid)
    del mid, res[0], content
    Ndata = []
    for i in range(0, len(res)):
        try:
            Ndataitem = {}
            k = 0

            if res[i][0][0] in ["R", "S"]:
            # print(res[i][0][0])
            # if res[i][0][0] == "R":
                for j in range(0, len(res[i])):
                    Ndataitem['数据组号'] = i + 1
                    if j == 0:
                        Ndataitem['卫星PRN号'] = res[i][j].split('\n')[0:3][0][0:3].replace(" ", '0')
                        Ndataitem['历元'] = res[i][j].strip("\n")[3:22].replace("  ", ' ')
                        if isin(Ndataitem['卫星PRN号'], Ndataitem['历元'], Ndata):
                            # 由于卫星广播星历更新时间为2个小时所以n文件内大部分数据冗余，此处判断若同一历元，同一卫星的星历数据已经存入，则不再存入
                            break
                        Ndataitem['svc'] = float(str(res[i][j].strip('\n')[23:42]).replace(" ", '').replace("D", "e"))
                        Ndataitem['svr'] = float(
                            str(res[i][j].strip("\n")[42:41 + 20]).replace(" ", '').replace("D", "e"))
                        Ndataitem['TEO'] = float(
                            str(res[i][j].strip("\n")[61:]).replace(" ", '').replace("D", "e"))
                    else:
                        try:
                            Ndataitem[glonassitemkeys[k]] = float(
                                str(res[i][j].strip('\n'))[3:23].replace(" ", '').replace("D", "e")) * 1000
                        except:
                            Ndataitem[glonassitemkeys[k]] = 0
                        try:
                            Ndataitem[glonassitemkeys[k + 1]] = float(
                                str(res[i][j].strip('\n'))[23:42].replace(" ", '').replace("D", "e")) * 1000
                        except:
                            Ndataitem[glonassitemkeys[k + 1]] = 0
                        try:
                            Ndataitem[glonassitemkeys[k + 2]] = float(
                                str(res[i][j].strip('\n'))[42:61].replace(" ", '').replace("D", "e")) * 1000
                        except:
                            Ndataitem[glonassitemkeys[k + 2]] = 0
                        try:
                            Ndataitem[glonassitemkeys[k + 3]] = float(
                                str(res[i][j].strip('\n'))[61:].replace(" ", '').replace("D", "e"))
                        except:
                            Ndataitem[glonassitemkeys[k + 3]] = 0
                        k = k + 4

                if len(Ndataitem) > 4:
                    Ndata.append(Ndataitem)
            else:
                try:
                    Ndataitem['数据组号'] = i + 1
                    for j in range(0, len(res[i])):
                        if j == 0:
                            Ndataitem['卫星PRN号'] = res[i][j].split('\n')[0:3][0][0:3].replace(" ", '')
                            Ndataitem['历元'] = res[i][j].strip("\n")[3:22].replace("  ", ' ')
                            if isin(Ndataitem['卫星PRN号'], Ndataitem['历元'], Ndata):
                                # 由于卫星广播星历更新时间为2个小时所以n文件内大部分数据冗余，此处判断若同一历元，同一卫星的星历数据已经存入，则不再存入
                                break
                            Ndataitem['卫星钟差（s）'] = float(
                                str(res[i][j].strip('\n')[23:42]).replace(" ", '').replace("D", "e"))
                            Ndataitem['卫星钟漂（s/s）'] = float(
                                str(res[i][j].strip("\n")[42:41 + 20]).replace(" ", '').replace("D", "e"))
                            Ndataitem['卫星钟漂速度参数（s/s/s）'] = float(
                                str(res[i][j].strip("\n")[61:]).replace(" ", '').replace("D", "e"))
                        else:
                            try:
                                Ndataitem[Ndataitemkeys[k]] = float(
                                    str(res[i][j].strip('\n'))[3:23].replace(" ", '').replace("D", "e"))
                            except:
                                print(k)
                                Ndataitem[Ndataitemkeys[k]] = 0
                            try:
                                Ndataitem[Ndataitemkeys[k + 1]] = float(
                                    str(res[i][j].strip('\n'))[23:42].replace(" ", '').replace("D", "e"))
                            except:
                                Ndataitem[Ndataitemkeys[k + 1]] = 0
                            try:
                                Ndataitem[Ndataitemkeys[k + 2]] = float(
                                    str(res[i][j].strip('\n'))[42:61].replace(" ", '').replace("D", "e"))
                            except:
                                Ndataitem[Ndataitemkeys[k + 2]] = 0
                            try:
                                Ndataitem[Ndataitemkeys[k + 3]] = float(
                                    str(res[i][j].strip('\n'))[61:].replace(" ", '').replace("D", "e"))
                            except:
                                Ndataitem[Ndataitemkeys[k + 3]] = 0
                            k = k + 4
                except:
                    pass
                if len(Ndataitem) > 4:
                    Ndata.append(Ndataitem)
        except:
            pass
    return Ndata


def space_num(odata, type, stype):
    # n = stype[type]['all']
    n = 32
    content = []
    mid = 0
    # print(odata)
    for i in range(0, n):
        if bt.all_spece(odata[0:12]):
            mid = 0
            odata = odata[16:]
        else:
            # print(odata[0:12].replace(' ',''))
            mid = float(odata[0:12])
            odata = odata[16:]
        content.append(mid)
    del odata, n, mid, stype
    return content


def gloposition(x, y, z, u, v, w, xls, yls, zls, h):
    # 四阶龙格库塔

    def glonassf(xi, yi, zi, vx, vy, vz, xls, yls, zls):
        GM = 3.986008e14
        r = math.sqrt(xi ** 2 + yi ** 2 + zi ** 2)
        ae = 6.378137e6
        C02 = 1.08263e-3
        we = 7.292115e-5
        ax = -GM * xi / (r ** 3) - 1.5 * C02 * GM * ae * ae * xi * (1 - (5 * (zi ** 2)) / (r ** 2)) / (
                r ** 5) + we * we * xi + 2 * we * vy + xls
        ay = -GM * yi / (r ** 3) - 1.5 * C02 * GM * ae * ae * yi * (1 - (5 * (zi ** 2)) / (r ** 2)) / (
                r ** 5) + we * we * yi - 2 * we * vx + yls
        az = -GM * zi / (r ** 3) - 1.5 * C02 * GM * ae * ae * zi * (3 - (5 * (zi ** 2)) / (r ** 2)) / (r ** 5) + zls

        return [ax, ay, az]

    def rk4(inti, step):
        x, y, z, u, v, w, xls, yls, zls = inti
        a1 = glonassf(x, y, z, u, v, w, xls, yls, zls)
        u1 = u
        v1 = v
        w1 = w
        a2 = glonassf(x + u1 * h / 2, y + v1 * h / 2, z + w1 * h / 2, u + a1[0] * h / 2, v + a1[1] * h / 2,
                      w + a1[2] * h / 2, xls, yls, zls)
        u2 = u + a1[0] * h / 2
        v2 = v + a1[1] * h / 2
        w2 = w + a1[2] * h / 2
        a3 = glonassf(x + u2 * h / 2, y + v2 * h / 2, z + w2 * h / 2, u + a2[0] * h / 2, v + a2[1] * h / 2,
                      w + a2[2] * h / 2, xls, yls, zls)
        u3 = u + a2[0] * h / 2
        v3 = v + a2[1] * h / 2
        w3 = w + a2[2] * h / 2
        a4 = glonassf(x + u3 * h, y + v3 * h, z + w3 * h, u + a3[0] * h, v + a3[1] * h, w + a3[2] * h, xls, yls, zls)
        u4 = u + a3[0] * h
        v4 = v + a3[1] * h
        w4 = w + a3[2] * h
        u = u + h * (a1[0] + 2 * a2[0] + 2 * a3[0] + a4[0]) / 6
        v = v + h * (a1[1] + 2 * a2[1] + 2 * a3[1] + a4[1]) / 6
        w = w + h * (a1[2] + 2 * a2[2] + 2 * a3[2] + a4[2]) / 6
        x = x + h * (u1 + 2 * u2 + 2 * u3 + u4) / 6
        y = y + h * (v1 + 2 * v2 + 2 * v3 + v4) / 6
        z = z + h * (w1 + 2 * w2 + 2 * w3 + w4) / 6

        return [x, y, z, u, v, w, a4[0], a4[1], a4[2]]
    #mm = mt.arithmetic_sequence(0,h,200)

    mid = [x, y, z, u, v, w, xls, yls, zls]

    mid = rk4(mid,1)
    return mid[0:3]




@bt.timeshow
def read_ofile(path:str):
    """
    docstrings
    RINEX VERSION / TYPE        :        rinex文件版本 / 文件类型: 本文件为O
    PGM / RUN BY / DATE         :        创建本数据文件所采用程序的名称 / 创建本数据文件单位的名称 / 创建本数据文件的日期
    MARKER NUMBER               :        天线标记名称
    OBSERVER / AGENCY           :        观测员姓名 / 观测单位名称
    REC # / TYPE / VERS         :        接收器编号、类型和版本
    ANT # / TYPE                :        天线序列号、天线类型
    APPROX POSITION XYZ         :        测站的近似位置X、Y、Z
    ANTENNA: DELTA H/E/N        :        天线高、天线中心相对于测站标志在东向偏移量、北向偏移量
    SYS / # / OBS TYPES         :        观测类型
    INTERVAL                    :        观测历元间隔（频率）
    TIME OF FIRST OBS           ：       首次观测时间
    TIME OF LAST OBS            ：       终止观测时间
    GLONASS SLOT / FRQ #        ：
    RCV CLOCK OFFS APPL         :
    COMMENT                     :        注释
    """

    #加载o文件
    f :list= open(path).readlines()

    #头文件标签，创建头文件标签字典，值为空列表，准备存放头文件各类值
    labels :dict = {"RINEX VERSION / TYPE": [], "PGM / RUN BY / DATE": [], "MARKER NAME": [],
              "MARKER NUMBER": [], "OBSERVER / AGENCY": [],
              "REC # / TYPE / VERS": [], "ANT # / TYPE": [],
              "APPROX POSITION XYZ": [], "ANTENNA: DELTA H/E/N": [],
              "SYS / # / OBS TYPES": [],
              "INTERVAL": [], "TIME OF FIRST OBS": [],
              "TIME OF LAST OBS": [], "GLONASS SLOT / FRQ #": [],
              "RCV CLOCK OFFS APPL": [],
              "LEAP SECONDS": [], "# OF SATELLITES": [], "COMMENT": [],
              "PRN / # OF OBS": [], "SYS / PHASE SHIFT": [], "GLONASS COD/PHS/BIS": [],
              "MARKER TYPE": [], "SIGNAL STRENGTH UNIT": []}
    #获取头文件结束行
    lines_of_head:int = head_num(f)

    for i in range(lines_of_head - 1):
        #根据头文件标签，填充label
        labels[f[i][60:].rstrip("\n").rstrip(" ")].append(f[i][:60].strip(" "))


    #s_type 准备存贮观测类型及数量
    s_type :dict= {"G": {"C": [], "L": [], "S": [], "D": [], "all": 0},
              "R": {"C": [], "L": [], "S": [], "D": [], "all": 0},
              "C": {"C": [], "L": [], "S": [], "D": [], "all": 0},
              "E": {"C": [], "L": [], "S": [], "D": [], "all": 0},
              "J": {"C": [], "L": [], "S": [], "D": [], "all": 0},
              "I": {"C": [], "L": [], "S": [], "D": [], "all": 0},
              "S": {"C": [], "L": [], "S": [], "D": [], "all": 0}}
    sys_res :list= []
    raw :int= 0

    for i in range(len(labels['SYS / # / OBS TYPES'])):
        de = bt.separate(labels['SYS / # / OBS TYPES'][i])
        if de[0] in ['G', 'R', 'E', 'J', 'C', 'I', 'S']:
            sys_res.append([de[0], de[2:]])
            raw = raw + 1
        else:
            sys_res[raw - 1][1] = sys_res[raw - 1][1] + de
    bt.list_read(sys_res)
    del  raw
    for i in range(len(sys_res)):
        for j in range(len(sys_res[i][1])):
            s_type[sys_res[i][0]][sys_res[i][1][j][0]].append(j)
    for key in s_type.keys():
        s_type[key]['all'] = len(s_type[key]['C']) + len(s_type[key]['L']) + len(s_type[key]['S']) + len(
            s_type[key]['D'])
    content = f[lines_of_head:]
    cal, sep = [], []
    for i in range(0, len(content)):
        if len(content[i]) < 260:
            content[i] = content[i] + (260 - len(content[i])) * " "
        if content[i][0] == '>':
            cal.append(sep)
            sep = []
        sep.append(content[i])
    cal.append(sep)
    del sep, cal[0], content
    for i in range(0, len(cal)):
        cal[i][0] = cal[i][0][1:-1]
        cal[i][0] = bt.all_float_list(cal[i][0][:-1])
        #all_float_list，将列表中所有元素转换成浮点数
        for j in range(1, len(cal[i])):
            cal[i][j] = cal[i][j][:-2]
    for i in range(0, len(cal)):
        cal[i][0] = bt.all_float_list(bt.separate(cal[i][0][1:-2])[:-1])
        #separate，将一段字符串，按照“ ”分开，不用strip的原因是，strip返回列表里有“ ”
        for j in range(1, len(cal[i])):
            res = [cal[i][j][:3]]
            t = s_type[cal[i][j][0][0]]['all']
            for z in range(0, t):
                mid = cal[i][j][z * 16 + 4:(z + 1) * 16 + 4]
                if bt.all_space(mid):
                    #all_space判断是否为空值、空格、tab
                    mid = "0"
                mid = bt.separate(mid)
                if len(mid) == 2:
                    del mid[-1]
                try:
                    mid = float(mid[0])
                except:
                    mid = 0
                res.append(mid)
            cal[i][j] = res
    return cal, s_type, labels


@bt.timeshow
def choice(odata, s_type, target):
    """
    odata为观测值
    s_type为各观测值类型所在位置
    target为筛选目标系统与信号类型
    """
    """
    单个odata数据，如下
    [[2023.0, 3.0, 3.0, 6.0, 41.0, 13.0, 0.0], 
    ['G01', 21038457.241, 21038462.393, 21038463.126, 110557838.229, 86148986.644, 82559447.517, 46.7, 45.83, 51.74, 1322.366, 1030.53, 987.44], 
    ['G07', 20690779.219, 20690781.411, 0.0, 108730778.938, 84725289.714, 0.0, 46.48, 45.95, 0.0, -2122.866, -1654.137, 0.0], 
    ['G08', 21652744.104, 21652750.255, 21652748.872, 113785937.065, 88664391.725, 84970036.568, 46.03, 46.31, 50.77, -3527.878, -2748.963, -2634.483],
    ['G10', 25623858.191, 25623869.471, 25623868.46, 134654287.096, 104925442.123, 100553565.799, 39.08, 39.4, 42.14, -1221.448, -951.659, -912.086], 
    ['G14', 23296445.83, 23296453.849, 23296459.139, 122423648.636, 95395076.424, 91420308.245, 37.45, 41.36, 46.13, 1467.103, 1143.276, 1095.501],
    ['G16', 25180593.618, 25180602.128, 0.0, 132324914.001, 103110348.374, 0.0, 37.4, 33.7, 0.0, -4166.514, -3246.671, 0.0], 
    ['G17', 24127217.547, 24127228.218, 0.0, 126789387.107, 98796957.355, 0.0, 41.35, 38.98, 0.0, 1239.702, 965.878, 0.0], 
    ['G21', 20509714.727, 20509715.72, 0.0, 107779277.418, 83983855.5, 0.0, 46.8, 44.98, 0.0, -1507.914, -1175.009, 0.0], 
    ['G27', 24186324.086, 24186330.401, 24186333.836, 127099989.102, 99038989.883, 94912364.669, 39.02, 37.93, 44.02, -4469.855, -3482.938, -3337.92], 
    ['G30', 21829098.293, 21829103.51, 21829104.659, 114712688.523, 89386536.461, 85662095.099, 45.33, 45.23, 50.08, 197.655, 154.063, 147.582], 
    ['R09', 20831206.851, 20831210.844, 111237481.462, 86518061.243, 50.19, 49.02, -3963.091, -3082.38], 
    ['R10', 19406951.319, 0.0, 103449969.646, 0.0, 47.73, 0.0, -791.276, 0.0], 
    ['R11', 22082097.764, 22082104.617, 118000039.955, 91777832.29, 45.38, 45.27, 2285.747, 1777.796], 
    ['R20', 19372715.11, 19372722.864, 103594610.728, 80573619.32, 39.88, 43.57, -2278.239, -1771.759], 
    ['R21', 21550385.165, 21550393.567, 115320465.901, 89693729.801, 46.78, 46.96, 1998.281, 1554.314], 
    ['C01', 37806889.744, 37806894.302, 37806894.805, 196870396.293, 152232697.316, 159973344.817, 42.9, 46.2, 41.77, -1261.511, -975.454, -1025.025, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C02', 38624128.239, 38624131.127, 38624132.271, 201125971.297, 155523371.009, 163431338.618, 40.81, 47.05, 42.9, -1292.291, -999.258, -1049.975, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C03', 37355787.145, 37355789.491, 37355789.13, 194521389.13, 150416284.009, 158064569.464, 45.23, 47.07, 43.44, -1263.625, -977.065, -1026.813, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C04', 38784632.074, 38784635.836, 38784634.702, 201961765.68, 156169660.002, 164110482.927, 40.69, 45.97, 41.61, -1266.996, -979.673, -1029.496, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C05', 40068526.377, 40068537.653, 40068535.233, 208647329.453, 161339400.022, 169543086.013, 40.68, 42.68, 37.17, -1297.357, -1003.207, -1054.315, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C06', 38528238.63, 38528243.42, 38528241.18, 200626641.294, 155137272.103, 163025595.849, 41.92, 44.77, 42.31, -3119.015, -2411.73, -2534.336, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C07', 36226733.279, 36226731.797, 36226730.045, 188642106.745, 145870037.067, 153287149.498, 47.94, 49.41, 47.33, -1178.548, -911.351, -957.588, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C09', 37174201.15, 37174206.108, 37174202.141, 193575815.575, 149685123.556, 157296214.305, 43.0, 48.2, 43.39, -2827.734, -2186.545, -2297.761, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C10', 36890666.396, 36890669.53, 36890666.334, 192099380.499, 148543438.481, 156096481.854, 46.09, 49.1, 46.13, -1056.936, -817.222, -858.784, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C16', 38572557.543, 38572567.427, 38572563.209, 200857426.412, 155315744.527, 163213138.922, 43.82, 43.96, 40.73, -3111.844, -2406.282, -2528.524, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C27', 25204970.149, 0.0, 25204975.696, 131248889.084, 0.0, 106650498.146, 43.82, 0.0, 44.63, -4223.228, 0.0, -3431.675, 25204970.075, 25204979.502, 132453005.07, 98909743.341, 38.17, 46.71, -4261.881, -3182.585], 
    ['C29', 22897263.388, 0.0, 22897266.36, 119232058.322, 0.0, 96885826.124, 49.41, 0.0, 49.37, 631.868, 0.0, 513.436, 22897263.253, 22897269.084, 120325931.101, 89853801.564, 43.64, 49.83, 637.621, 476.116], 
    ['C30', 21899667.595, 0.0, 21899672.477, 114037315.512, 0.0, 92664679.309, 51.48, 0.0, 50.28, -2337.645, 0.0, -1899.481, 21899667.705, 21899674.096, 115083531.11, 85939024.086, 46.3, 50.62, -2359.068, -1761.657], 
    ['C36', 21876526.757, 0.0, 21876534.377, 113916814.256, 0.0, 92566774.935, 52.63, 0.0, 50.8, -1679.339, 0.0, -1364.548, 21876526.942, 21876537.696, 114961923.298, 85848231.412, 47.28, 51.06, -1694.719, -1265.56], 
    ['C38', 39776650.785, 0.0, 39776658.193, 207127461.374, 0.0, 168308059.214, 43.4, 0.0, 42.98, 614.424, 0.0, 499.233, 39776651.649, 39776664.06, 209027714.959, 156092172.667, 38.91, 42.78, 619.972, 462.984], 
    ['C39', 39613705.835, 0.0, 39613711.633, 206278959.418, 0.0, 167618577.482, 43.37, 0.0, 41.59, -3165.041, 0.0, -2571.747, 39613705.758, 39613717.666, 208171432.6, 155452734.312, 38.81, 43.78, -3193.98, -2385.122], 
    ['C40', 36364610.432, 0.0, 36364611.796, 189360068.11, 0.0, 153870573.86, 50.08, 0.0, 49.77, -1321.814, 0.0, -1074.097, 36364610.193, 36364614.252, 191097317.104, 142702556.284, 45.3, 50.44, -1333.933, -996.156], 
    ['C45', 23722468.662, 0.0, 23722466.302, 123529115.051, 0.0, 100377517.244, 46.16, 0.0, 47.25, 1129.865, 0.0, 918.161, 23722470.343, 23722469.607, 124662420.423, 93092063.177, 42.44, 47.97, 1140.227, 851.494], 
    ['C46', 24098173.472, 0.0, 24098172.622, 125485515.682, 0.0, 101967253.988, 44.76, 0.0, 46.54, -4038.463, 0.0, -3281.555, 24098175.063, 24098174.808, 126636765.617, 94566410.232, 41.34, 47.83, -4075.496, -3043.38], 
    ['C56', 36196850.858, 0.0, 36196854.671, 188486499.08, 0.0, 153160738.501, 48.52, 0.0, 46.81, -961.42, 0.0, -781.186, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C59', 37893827.116, 0.0, 37893818.634, 197323102.693, 0.0, 160341147.465, 46.04, 0.0, 46.12, -1300.626, 0.0, -1056.867, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['C60', 38568611.47, 0.0, 38568602.4, 200836877.82, 0.0, 163196376.771, 45.11, 0.0, 46.27, -1247.112, 0.0, -1013.384, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
    ['E14', 20145591.817, 20145596.367, 20145594.153, 105865800.712, 79055645.36, 81117957.149, 47.44, 51.26, 52.0, 2190.834, 1636.031, 1678.685], 
    ['E15', 25744471.665, 25744480.313, 25744477.436, 135288102.515, 101026869.595, 103662342.682, 40.13, 47.36, 46.85, -3363.593, -2511.726, -2577.235], 
    ['E19', 29075565.503, 29075579.958, 29075576.035, 152793139.459, 114098819.043, 117075297.437, 29.33, 35.88, 33.93, -3601.173, -2689.264, -2759.43], 
    ['E21', 23860957.908, 23860963.309, 23860960.628, 125390178.828, 93635544.454, 96078201.379, 45.53, 50.74, 50.07, -457.186, -341.427, -350.313], 
    ['E26', 26060223.209, 26060233.39, 26060230.999, 136947397.851, 102265954.004, 104933748.998, 40.21, 44.89, 45.92, 1605.572, 1199.055, 1230.272], 
    ['E27', 24961582.07, 24961589.309, 24961587.343, 131174000.827, 97954637.555, 100509967.502, 43.74, 49.13, 49.2, -2644.192, -1974.604, -2026.069], 
    ['J02', 38608684.091, 38608688.436, 38608689.698, 202890006.062, 158096123.697, 151508791.173, 45.95, 47.77, 50.16, -592.154, -461.392, -442.187], 
    ['J03', 38503806.383, 38503810.672, 38503812.338, 202338866.432, 157666662.022, 151097230.548, 46.13, 48.04, 51.07, -1841.895, -1435.183, -1375.419]]
  
    
    """
    """
    s_type
    {
    'G': {'C': [0, 1, 2], 'L': [3, 4, 5], 'S': [6, 7, 8], 'D': [9, 10, 11], 'all': 12}, 
    'R': {'C': [0, 1], 'L': [2, 3], 'S': [4, 5], 'D': [6, 7], 'all': 8}, 
    'C': {'C': [0, 1, 2, 12, 13], 'L': [3, 4, 5, 14, 15], 'S': [6, 7, 8, 16, 17], 'D': [9, 10, 11, 18, 19], 'all': 20}, 
    'E': {'C': [0, 1, 2], 'L': [3, 4, 5], 'S': [6, 7, 8], 'D': [9, 10, 11], 'all': 12}, 
    'J': {'C': [0, 1, 2], 'L': [3, 4, 5], 'S': [6, 7, 8], 'D': [9, 10, 11], 'all': 12}, 
    'I': {'C': [], 'L': [], 'S': [], 'D': [], 'all': 0}, 
    'S': {'C': [], 'L': [], 'S': [], 'D': [], 'all': 0}
    }
    """
    mid, res, small = [], [], []
    #逐历元
    for i in range(0, len(odata)):
        #时间，年月日分秒
        time = odata[i][0][:-1]
        name = []
        mid = []
        #逐卫星
        for j in range(1, len(odata[i])):
            small = []
            """
            odata[i][j][0][0] 为卫星类型 "G C R E J"
            target[0]为目标卫星类型，可以多个
            """
            if odata[i][j][0][0] in target[0]:
                name.append(odata[i][j][0])
                
                """
                target[1]为目标信号类型
                """
                for x in range(0, len(target[1])):
                    for y in s_type[odata[i][j][0][0]][target[1][x]]:
                        small.append(odata[i][j][y + 1])
                mid.append(small)
        res.append([time, name, mid])
    """
    单个res如下所示：target为(["C"],["S"])
    [[2023.0, 3.0, 3.0, 6.0, 41.0, 13.0], 
    ['C01', 'C02', 'C03', 'C04', 'C05', 'C06', 'C07', 'C09', 'C10', 'C16', 'C27', 'C29', 'C30', 'C36', 'C38', 'C39', 'C40', 'C45', 'C46', 'C56', 'C59', 'C60'], 
    [[42.9, 46.2, 41.77, 0.0, 0.0], 
    [40.81, 47.05, 42.9, 0.0, 0.0], 
    [45.23, 47.07, 43.44, 0.0, 0.0], 
    [40.69, 45.97, 41.61, 0.0, 0.0], 
    [40.68, 42.68, 37.17, 0.0, 0.0], 
    [41.92, 44.77, 42.31, 0.0, 0.0], 
    [47.94, 49.41, 47.33, 0.0, 0.0], 
    [43.0, 48.2, 43.39, 0.0, 0.0], 
    [46.09, 49.1, 46.13, 0.0, 0.0], 
    [43.82, 43.96, 40.73, 0.0, 0.0], 
    [43.82, 0.0, 44.63, 38.17, 46.71], 
    [49.41, 0.0, 49.37, 43.64, 49.83], 
    [51.48, 0.0, 50.28, 46.3, 50.62], 
    [52.63, 0.0, 50.8, 47.28, 51.06], 
    [43.4, 0.0, 42.98, 38.91, 42.78], 
    [43.37, 0.0, 41.59, 38.81, 43.78], 
    [50.08, 0.0, 49.77, 45.3, 50.44], 
    [46.16, 0.0, 47.25, 42.44, 47.97], 
    [44.76, 0.0, 46.54, 41.34, 47.83], 
    [48.52, 0.0, 46.81, 0.0, 0.0], 
    [46.04, 0.0, 46.12, 0.0, 0.0], 
    [45.11, 0.0, 46.27, 0.0, 0.0]]]
    """
    return res


@bt.timeshow
def choicek(odata, s_type, target):
    mid, res, small = [], [], []
    for i in range(0, len(odata)):
        time = odata[i][0][:]
        name = []
        mid = []
        for j in range(1, len(odata[i])):
            small = []
            if odata[i][j][0][0] in target[0]:
                name.append(odata[i][j][0])
                for x in range(0, len(target[1])):
                    for y in s_type[odata[i][j][0][0]][target[1][x]]:
                        small.append(odata[i][j][y + 1])
                mid.append(small)
        res.append([time, name, mid])
    return res


def julday(cal):
    # 计算儒略日
    mjd = cal
    y = mjd[0]
    m = mjd[1]
    d = mjd[2]

    try:
        h = (mjd[3] * 3600 + mjd[4] * 60 + mjd[5]) / 3600
    except:
        h = (mjd[3] * 3600 + mjd[4] ) / 3600
    if m <= 2:
        y = y - 1
        m = m + 12
    jd = floor(365.25 * (y + 4716)) + floor(30.6001 * (m + 1)) + d + (h / 24) - 1537.5
    return jd


def gpstime(julday):
    a = floor(julday + 0.5)
    b = a + 1537
    c = floor((b - 122.1) / 365.25)
    e = floor(365.25 * c)
    f = floor((b - e) / 30.6001)
    d = b - e - floor(30.6001 * f) + (julday + 0.5) % 1
    day_of_week = floor(julday + 0.5) % 7
    sec_of_week = (d % 1 + day_of_week + 1) * 86400
    return sec_of_week


def find_eph(eph, sats, time):
    """
    eph为提取的n文件数据，sats为历元内单个卫星号，time为周秒
    """
    if sats[0] == "R":
        isat = []
        for i in range(0, len(eph)):
            if eph[i]['卫星PRN号'] == sats:
                isat.append(i)
        n = len(isat)
        if n == 0:
            icol = 0
            return 0
        icol = isat[0]
        glotime1 = bt.all_float_list(bt.separate(eph[icol]['历元']))
        glotime1 = gpstime(julday(glotime1))
        dtmin = glotime1 - time
        for t in isat:
            glotime = bt.all_float_list(bt.separate(eph[t]['历元']))
            glotime = gpstime(julday(glotime))
            dt = glotime - time
            # print(dt/60)
            if abs(dt) < abs(dtmin):
                icol = t
                dtmin = dt
        return icol
    else:
        isat = []
        for i in range(0, len(eph)):
            if eph[i]['卫星PRN号'] == sats:
                isat.append(i)
        n = len(isat)
        if n == 0:
            icol = 0
            return 0
        icol = isat[0]
        dtmin = eph[icol]['TEO'] - time
        for t in isat:
            dt = eph[t]['TEO'] - time
            if dt < 0:
                if abs(dt) < abs(dtmin):
                    icol = t
                    dtmin = dt
        return icol


def all_ofile_satname(o_choice_res):
    """
    输入选择后的o文件数据，返回每个历元出现的卫星号,以及每个历元的周秒
    """
    satlin = []
    sow = []
    for i in range(0, len(o_choice_res)):
        satlin.append(o_choice_res[i][1])
        jd = julday(o_choice_res[i][0])
        sow.append(gpstime(jd))
    return satlin, sow


def check_t(time):
    """
    修复GPS时间的上溢和下溢
    """
    half_week = 302400
    tt = time
    if time > half_week:
        tt = time - 2 * half_week
    if time < -half_week:
        tt = time + 2 * half_week
    return tt


def satpos(t, eph):
    satp = []
    GM = 3.986008e14
    Omegae_dot = 7.2921151467e-5
    roota = eph['sqrt_A']
    toe = eph['TEO']
    deltan = eph['Delta_n']
    M0 = eph['M0']
    i0 = eph['i0']
    omega = eph['omega']
    ecc = eph['e']
    cuc = eph['Cuc']
    cus = eph['Cus']
    crc = eph['Crc']
    crs = eph['Crs']
    cic = eph['Cic']
    cis = eph['Cis']
    omega0 = eph['OMEGA_A0']
    omegadot = eph['OMEGA_DOT']
    idot = eph['IDOT']
    A = roota * roota
    tk = check_t(t - toe)
    n0 = sqrt(GM / (A ** 3))
    n = n0 + deltan
    M = M0 + n * tk
    M = (M + 2 * pi) % (2 * pi)
    E = M
    for i in range(0, 10):
        E_old = E
        E = M + eph['e'] * sin(E)
        dE = (E - E_old) % (2 * pi)
        if abs(dE) < 1e-12:
            break
    E = (E + 2 * pi) % (2 * pi)
    v = atan2(sqrt(1 - ecc ** 2) * sin(E), cos(E) - ecc)
    phi = v + omega
    phi = (phi) % (2 * pi)
    u = phi + cuc * cos(2 * phi) + cus * sin(2 * phi)
    r = A * (1 - ecc * cos(E)) + crc * cos(2 * phi) + crs * sin(2 * phi)
    i = i0 + idot * tk + cic * cos(2 * phi) + cis * sin(2 * phi)
    omega = omega0 + (omegadot - Omegae_dot) * tk - Omegae_dot * toe
    omega = (omega + 2 * pi) % (2 * pi)
    x1 = cos(u) * r
    y1 = sin(u) * r
    satp.append(x1 * cos(omega) - y1 * cos(i) * sin(omega))
    satp.append(x1 * sin(omega) + y1 * cos(i) * cos(omega))
    satp.append(y1 * sin(i))
    return satp


def satposgeo(t, eph):
    satp = []
    GM = 3.986008e14
    Omegae_dot = 7.2921151467e-5
    roota = eph['sqrt_A']
    toe = eph['TEO']
    deltan = eph['Delta_n']
    M0 = eph['M0']
    i0 = eph['i0']
    omega = eph['omega']
    ecc = eph['e']
    cuc = eph['Cuc']
    cus = eph['Cus']
    crc = eph['Crc']
    crs = eph['Crs']
    cic = eph['Cic']
    cis = eph['Cis']
    omega0 = eph['OMEGA_A0']
    omegadot = eph['OMEGA_DOT']
    idot = eph['IDOT']
    A = roota * roota
    tk = check_t(t - toe)
    n0 = sqrt(GM / (A ** 3))
    n = n0 + deltan
    M = M0 + n * tk
    M = (M + 2 * pi) % (2 * pi)
    E = M
    for i in range(0, 10):
        E_old = E
        E = M + eph['e'] * sin(E)
        dE = (E - E_old) % (2 * pi)
        if abs(dE) < 1e-12:
            break
    E = (E + 2 * pi) % (2 * pi)
    v = atan2(sqrt(1 - ecc ** 2) * sin(E), cos(E) - ecc)
    phi = v + omega
    phi = (phi) % (2 * pi)
    u = phi + cuc * cos(2 * phi) + cus * sin(2 * phi)
    r = A * (1 - ecc * cos(E)) + crc * cos(2 * phi) + crs * sin(2 * phi)
    i = i0 + idot * tk + cic * cos(2 * phi) + cis * sin(2 * phi)
    omega = omega0 + (omegadot - Omegae_dot) * tk - Omegae_dot * toe
    omega = (omega + 2 * pi) % (2 * pi)
    x1 = cos(u) * r
    y1 = sin(u) * r
    xgk = x1 * cos(omega) - y1 * cos(i) * sin(omega)
    ygk = x1 * sin(omega) + y1 * cos(i) * cos(omega)
    zgk = y1 * sin(i)
    satp.append(cos(omegadot * tk) * xgk + sin(omegadot * tk) * cos(-5 / 180 * pi) * ygk + sin(omegadot * tk) * sin(
        -5 / 180 * pi) * zgk)
    satp.append(-sin(omegadot * tk) * xgk + cos(omegadot * tk) * cos(-5 / 180 * pi) * ygk + cos(omegadot * tk) * sin(
        -5 / 180 * pi) * zgk)
    satp.append(-sin(-5 / 180 * pi) * ygk + cos(-5 / 180 * pi) * zgk)
    return satp


def recpo_lsa(obs, sats, time, Eph):
    col_Eph = []
    El = []
    v_light = 299792458
    dtr = pi / 180

    m = len(obs)
    el = zeros((m, 1))
    for t in range(0, m):
        col = find_eph(Eph, sats[t], time)
        col_Eph.append(col)
    OBS = []
    for i in range(0, len(obs)):
        OBS.append(obs[i])
    xpos = []
    n = len(col_Eph)
    if n == 0:
        return nan
    for i in range(0, n):
        k = col_Eph[i]
        if sats[i][0] == "R":
            a = check_t(time)
            if a < 0:
                a = 86400 + a
            #h = check_t(a - Eph[k]['TEO']) + 3 * 3600 - 15 * 60
            h = check_t(a - Eph[k]['TEO']) + (1.5 * 3600+15*60+46)
            X = gloposition(Eph[k]['x'], Eph[k]['y'], Eph[k]['z'], Eph[k]['xv'], Eph[k]['yv'], Eph[k]['zv'],
                            Eph[k]['xa'], Eph[k]['ya'], Eph[k]['za'], h,)
            xpos.append(X[:3])
        else:
            try:
                tx_RAW = time - OBS[i] / v_light
            except:
                tx_RAW = time
            toc = Eph[k]['TEO']
            dt = check_t(tx_RAW - toc)
            try:
                tcorr = (Eph[k]['卫星钟漂速度参数（s/s/s）'] * dt + Eph[k]['卫星钟漂（s/s）']) * dt + Eph[k]['卫星钟差（s）']
            except:
                tcorr = 0
            tx_GPS = tx_RAW - tcorr
            dt = check_t(tx_GPS - toc)
            try:
                tcorr = (Eph[k]['卫星钟漂速度参数（s/s/s）'] * dt + Eph[k]['卫星钟漂（s/s）']) * dt + Eph[k]['卫星钟差（s）']
            except:
                tcorr = 0
            tx_GPS = tx_RAW - tcorr
            if sats[i] in ["C01", "C02", "C03", "C04", "C05", "C59", "C60"]:
                X = satposgeo(tx_GPS, Eph[k])
            else:
                X = satpos(tx_GPS, Eph[k])
            xpos.append(X)
    return xpos


def e_r_corr(traveltime, X_sat):
    # 返回信号传播时间内由于地球自转而旋转的卫星ECEF坐标
    # ECEF地心地固坐标系
    Omegae_dot = 7.292115147e-5
    omegatau = Omegae_dot * traveltime
    R3 = array(
        [[cos(omegatau), sin(omegatau), 0],
         [-sin(omegatau), cos(omegatau), 0],
         [0, 0, 1]])
    X_sat_rot = matmul(R3, X_sat)
    for i in range(0, len(X_sat_rot)):
        if type(X_sat_rot[i]) == 'numpy.ndarray':
            X_sat_rot[i] = X_sat_rot[i].tolist()[0]
    return X_sat_rot


def toBLH(x, y, z):
    L = atan(y / x)
    a = 6378137
    c = 6399593.6258
    e = 1 / 298.257223563
    e12 = 0.00673949674227
    dtr = 180 / pi
    t0 = z / (sqrt(x ** 2 + y ** 2))
    p = c * (e ** 2) / (sqrt(x ** 2 + y ** 2))
    k = 1 + e12
    t1 = t0 + p * t0 / (sqrt(k + t0 ** 2))
    t2 = t0
    while abs(t2 - t1) > 1e-8:
        t1 = t2
        t2 = t0 + p * t1 / (k + t1 ** 2)
    b = atan(t2)
    h = (z / sin(b)) - a / (sqrt(1 - (e * sin(b)) ** 2)) * (1 - e ** 2)
    return b, L + pi, h

def Blhtoxyz84(lat,lon,h):
    a = 6378137
    b = 6356752.31424518
    E = (a*a-b*b)/(a*a)
    dtr = math.pi/180
    coslat = math.cos(lat*dtr)
    sinlat = math.sin(lat*dtr)
    coslon = math.cos(lon*dtr)
    sinlon = math.sin(lon*dtr)
    N = a/(math.sqrt(1-E*sinlat*sinlat))
    NH = N+h
    x = NH*coslat*coslon
    y = NH*coslat*sinlon
    z = (b*b*N/(a**2)+h)*sinlat
    return [x,y,z]


@bt.timeshow
def az(opath, npath, target):
    v_light = 299792458
    dtr = pi / 180
    obs, s_type, label = read_ofile(opath)
    obs = choice(obs, s_type, target)
    satlin, sow = all_ofile_satname(obs)
    time = []
    for i in range(0, len(obs)):
        time.append(obs[i][0])
        try:
            del obs[i][0], obs[i][0]
        except:
            pass
    for i in range(0, len(obs)):
        obs[i] = obs[i][0]
    for i in range(0, len(obs)):
        for j in range(0, len(obs[i])):
            for z in range(1, len(obs[i][j])):
                del obs[i][j][-1]
    Eph = readnfile(npath)
    pos = bt.all_float_list(bt.separate(str(label['APPROX POSITION XYZ'][0])))
    en = len(obs)
    B, L, H = toBLH(pos[0], pos[1], pos[2])
    Transformationmatrix = array([[-sin(B) * cos(L), -sin(B) * sin(L), cos(B)],
                                  [-sin(L), cos(L), 0],
                                  [cos(B) * cos(L), cos(B) * sin(L), sin(B)]
                                  ])
    Surveystation = []
    mid = []
    for q in range(0, en):
        weeksow = sow[q]
        sats = satlin[q]
        a = recpo_lsa(obs[q], sats, weeksow, Eph)
        # a为卫星位置
        mid = []
        to = (180 / pi)
        for i in range(0, len(a)):
            cha = array([[a[i][0] - pos[0]], [a[i][1] - pos[1]], [a[i][2] - pos[2]]])
            tar = dot(Transformationmatrix, cha)
            tar = tar.tolist()
            tar = [tar[0][0], tar[1][0], tar[2][0]]
            # 地球自转改正
            traceltime = sqrt(tar[0] ** 2 + tar[1] ** 2 + tar[2] ** 2) / v_light
            sateloc = e_r_corr(traceltime, a[i])
            cha = array([[sateloc[0] - pos[0]], [sateloc[1] - pos[1]], [sateloc[2] - pos[2]]])
            tar = dot(Transformationmatrix, cha)
            tar = tar.tolist()
            tar = [tar[0][0], tar[1][0], tar[2][0]]
            r = sqrt(tar[0] ** 2 + tar[1] ** 2 + tar[2] ** 2)
            ele = 90 - arcsin(tar[2] / r) * to
            azi = arctan2(tar[1], tar[0])
            if ele > 90:
                mid.append([0, 0])
            else:
                mid.append([azi, ele, ])
        Surveystation.append(mid)
    return time, satlin, Surveystation


def diffentsatellite(allsatellite, obssatellite):
    x = []
    for i in allsatellite:
        if i not in obssatellite:
            x.append(i)
    return x

@bt.timeshow
def azinitialization(time, ae, satlin, type):
    satellite = {
        "G": {
            "01": [[], [], []], "02": [[], [], []], "03": [[], [], []], "04": [[], [], []], "05": [[], [], []],
            "06": [[], [], []], "07": [[], [], []], "08": [[], [], []], "09": [[], [], []], "10": [[], [], []],
            "11": [[], [], []], "12": [[], [], []], "13": [[], [], []], "14": [[], [], []], "15": [[], [], []],
            "16": [[], [], []], "17": [[], [], []], "18": [[], [], []], "19": [[], [], []], "20": [[], [], []],
            "21": [[], [], []], "22": [[], [], []], "23": [[], [], []], "24": [[], [], []], "25": [[], [], []],
            "26": [[], [], []], "27": [[], [], []], "28": [[], [], []], "29": [[], [], []], "30": [[], [], []],
            "31": [[], [], []], "32": [[], [], []], "33": [[], [], []], "34": [[], [], []], "35": [[], [], []],
            "36": [[], [], []], "37": [[], [], []], "38": [[], [], []], "39": [[], [], []], "40": [[], [], []],
            "41": [[], [], []], "42": [[], [], []], "43": [[], [], []], "44": [[], [], []], "45": [[], [], []],
            "46": [[], [], []], "47": [[], [], []], "48": [[], [], []], "49": [[], [], []], "50": [[], [], []],
            "51": [[], [], []], "52": [[], [], []], "53": [[], [], []], "54": [[], [], []], "55": [[], [], []],
            "56": [[], [], []], "57": [[], [], []], "58": [[], [], []], "59": [[], [], []], "60": [[], [], []], },
        "E": {
            "01": [[], [], []], "02": [[], [], []], "03": [[], [], []], "04": [[], [], []], "05": [[], [], []],
            "06": [[], [], []], "07": [[], [], []], "08": [[], [], []], "09": [[], [], []], "10": [[], [], []],
            "11": [[], [], []], "12": [[], [], []], "13": [[], [], []], "14": [[], [], []], "15": [[], [], []],
            "16": [[], [], []], "17": [[], [], []], "18": [[], [], []], "19": [[], [], []], "20": [[], [], []],
            "21": [[], [], []], "22": [[], [], []], "23": [[], [], []], "24": [[], [], []], "25": [[], [], []],
            "26": [[], [], []], "27": [[], [], []], "28": [[], [], []], "29": [[], [], []], "30": [[], [], []],
            "31": [[], [], []], "32": [[], [], []], "33": [[], [], []], "34": [[], [], []], "35": [[], [], []],
            "36": [[], [], []], "37": [[], [], []], "38": [[], [], []], "39": [[], [], []], "40": [[], [], []],
            "41": [[], [], []], "42": [[], [], []], "43": [[], [], []], "44": [[], [], []], "45": [[], [], []],
            "46": [[], [], []], "47": [[], [], []], "48": [[], [], []], "49": [[], [], []], "50": [[], [], []],
            "51": [[], [], []], "52": [[], [], []], "53": [[], [], []], "54": [[], [], []], "55": [[], [], []],
            "56": [[], [], []], "57": [[], [], []], "58": [[], [], []], "59": [[], [], []], "60": [[], [], []], },
        "C": {
            "01": [[], [], []], "02": [[], [], []], "03": [[], [], []], "04": [[], [], []], "05": [[], [], []],
            "06": [[], [], []], "07": [[], [], []], "08": [[], [], []], "09": [[], [], []], "10": [[], [], []],
            "11": [[], [], []], "12": [[], [], []], "13": [[], [], []], "14": [[], [], []], "15": [[], [], []],
            "16": [[], [], []], "17": [[], [], []], "18": [[], [], []], "19": [[], [], []], "20": [[], [], []],
            "21": [[], [], []], "22": [[], [], []], "23": [[], [], []], "24": [[], [], []], "25": [[], [], []],
            "26": [[], [], []], "27": [[], [], []], "28": [[], [], []], "29": [[], [], []], "30": [[], [], []],
            "31": [[], [], []], "32": [[], [], []], "33": [[], [], []], "34": [[], [], []], "35": [[], [], []],
            "36": [[], [], []], "37": [[], [], []], "38": [[], [], []], "39": [[], [], []], "40": [[], [], []],
            "41": [[], [], []], "42": [[], [], []], "43": [[], [], []], "44": [[], [], []], "45": [[], [], []],
            "46": [[], [], []], "47": [[], [], []], "48": [[], [], []], "49": [[], [], []], "50": [[], [], []],
            "51": [[], [], []], "52": [[], [], []], "53": [[], [], []], "54": [[], [], []], "55": [[], [], []],
            "56": [[], [], []], "57": [[], [], []], "58": [[], [], []], "59": [[], [], []], "60": [[], [], []],
        },
        "R": {
            "01": [[], [], []], "02": [[], [], []], "03": [[], [], []], "04": [[], [], []], "05": [[], [], []],
            "06": [[], [], []], "07": [[], [], []], "08": [[], [], []], "09": [[], [], []], "10": [[], [], []],
            "11": [[], [], []], "12": [[], [], []], "13": [[], [], []], "14": [[], [], []], "15": [[], [], []],
            "16": [[], [], []], "17": [[], [], []], "18": [[], [], []], "19": [[], [], []], "20": [[], [], []],
            "21": [[], [], []], "22": [[], [], []], "23": [[], [], []], "24": [[], [], []], "25": [[], [], []],
            "26": [[], [], []], "27": [[], [], []], "28": [[], [], []], "29": [[], [], []], "30": [[], [], []],
            "31": [[], [], []], "32": [[], [], []], "33": [[], [], []], "34": [[], [], []], "35": [[], [], []],
            "36": [[], [], []], "37": [[], [], []], "38": [[], [], []], "39": [[], [], []], "40": [[], [], []],
            "41": [[], [], []], "42": [[], [], []], "43": [[], [], []], "44": [[], [], []], "45": [[], [], []],
            "46": [[], [], []], "47": [[], [], []], "48": [[], [], []], "49": [[], [], []], "50": [[], [], []],
            "51": [[], [], []], "52": [[], [], []], "53": [[], [], []], "54": [[], [], []], "55": [[], [], []],
            "56": [[], [], []], "57": [[], [], []], "58": [[], [], []], "59": [[], [], []], "60": [[], [], []],
        },
        "J": {
            "01": [[], [], []], "02": [[], [], []], "03": [[], [], []], "04": [[], [], []], "05": [[], [], []],
            "06": [[], [], []], "07": [[], [], []], "08": [[], [], []], "09": [[], [], []], "10": [[], [], []],
            "11": [[], [], []], "12": [[], [], []], "13": [[], [], []], "14": [[], [], []], "15": [[], [], []],
            "16": [[], [], []], "17": [[], [], []], "18": [[], [], []], "19": [[], [], []], "20": [[], [], []],
            "21": [[], [], []], "22": [[], [], []], "23": [[], [], []], "24": [[], [], []], "25": [[], [], []],
            "26": [[], [], []], "27": [[], [], []], "28": [[], [], []], "29": [[], [], []], "30": [[], [], []],
            "31": [[], [], []], "32": [[], [], []], "33": [[], [], []], "34": [[], [], []], "35": [[], [], []],
            "36": [[], [], []], "37": [[], [], []], "38": [[], [], []], "39": [[], [], []], "40": [[], [], []],
            "41": [[], [], []], "42": [[], [], []], "43": [[], [], []], "44": [[], [], []], "45": [[], [], []],
            "46": [[], [], []], "47": [[], [], []], "48": [[], [], []], "49": [[], [], []], "50": [[], [], []],
            "51": [[], [], []], "52": [[], [], []], "53": [[], [], []], "54": [[], [], []], "55": [[], [], []],
            "56": [[], [], []], "57": [[], [], []], "58": [[], [], []], "59": [[], [], []], "60": [[], [], []],
        },
    }

    if type == 1:  # 笛卡尔
        for i in range(0, len(ae)):
            for j in range(0, len(ae[i])):
                m = ae[i][j][0]
                if m < 0:
                    m = m + 2 * math.pi
                satellite[satlin[i][j][0]][satlin[i][j][1:]][0].append(time[i])
                satellite[satlin[i][j][0]][satlin[i][j][1:]][1].append(m)
                satellite[satlin[i][j][0]][satlin[i][j][1:]][2].append(90 - ae[i][j][1])

    elif type == -1:  # 极坐标
        en =len(ae)
        for i in range(0, len(ae)):
            for j in range(0, len(ae[i])):
                satellite[satlin[i][j][0]][satlin[i][j][1:]][0].append(time[i])
                satellite[satlin[i][j][0]][satlin[i][j][1:]][1].append(ae[i][j][0])
                satellite[satlin[i][j][0]][satlin[i][j][1:]][2].append(ae[i][j][1])

    return satellite


@bt.timeshow
def drawaz(satellite, type):
    color = ["#660000", "#660066", "#6600CC", "#663300", "#663366", "#6633CC", "#666600", "#666666",
             "#6666CC", "#669900", "#669966", "#FF0000", "#FF0099", "#FF00FF", "#FF3366", "#FF33FF",
             "#FF6633", "#FF9900", "#FF9999", "#FFCC00", "#FFCC99", "#CC0066", "#33CCFF", "#996633",
             "#336666", "#33FFFF", "#669900", "#FF0033", "#000066", "#0000CC", "#9900FF", "#00FFCC",
             "#660000", "#660066", "#6600CC", "#663300", "#663366", "#6633CC", "#666600", "#666666",
             "#6666CC", "#669900", "#669966", "#FF0000", "#FF0099", "#FF00FF", "#FF3366", "#FF33FF",
             "#FF6633", "#FF9900", "#FF9999", "#FFCC00", "#FFCC99", "#CC0066", "#33CCFF", "#996633",
             "#336666", "#33FFFF", "#669900", "#FF0033", "#000066", "#0000CC", "#9900FF", "#00FFCC",
             ]
    # 笛卡尔坐标系
    if type == 1:
        for key in satellite.keys():
            for keys in satellite[key].keys():
                if len(satellite[key][keys][0]) > 1 and (key != "C" and keys not in ["56", "57", "58", "59"]):
                    plt.scatter(satellite[key][keys][1], satellite[key][keys][2], c=color[int(float(keys)) - 1],
                                marker=".",
                                alpha=0.75, s=5)
                    plt.annotate(key + keys, (satellite[key][keys][1][0], satellite[key][keys][2][0]))
        plt.xlabel("azimuth")
        plt.ylabel("elevation")
        plt.xlim(0, math.pi * 2)
        plt.ylim(0, 90)
        xticks = [0, 1 * math.pi / 3, 2 * math.pi / 3, math.pi, 4 * math.pi / 3, 5 * math.pi / 3, 2 * math.pi]
        xticklabel = ["0", "60", "120", "180", "240", "300", "360"]
        plt.xticks(xticks, xticklabel)
        plt.legend()
        plt.savefig("GNSS sky map.png", dpi=1000)
        plt.show()
    # 极坐标系
    elif type == -1:
        ax = plt.subplot(111, projection='polar')
        ax.set_theta_direction(-1)
        ax.set_theta_zero_location('N')
        ax.yaxis.set_label_position('right')
        ax.tick_params('y', labelleft=False)
        ax.grid(linestyle='--')
        labels = ['N', '45°', 'E', '135°', 'S', '225°', 'W', '315°']
        ax.set_thetagrids(range(0, 360, 45), labels, fontweight='semibold')
        for key in satellite.keys():
            for keys in satellite[key].keys():
                if len(satellite[key][keys][0]) > 1:
                    c = ax.scatter(satellite[key][keys][1], satellite[key][keys][2], c=color[int(float(keys)) - 1],
                                   marker=".",
                                   alpha=0.75, s=0.5)
                    ax.annotate(key + keys, (satellite[key][keys][1][0], satellite[key][keys][2][0]))

        SATAZ = [0, 0, 0, 0, 0, 0]
        SATEL = [0, 15, 30, 45, 60, 75]
        SATID = ['90', '75', '60', '45', '30', '15']
        for i in range(0, 6):
            ax.text(SATAZ[i], SATEL[i], SATID[i])
        ax.set_rticks(range(0, 105, 15))
        plt.legend()
        plt.savefig("laosang.png",dpi =1000)
        plt.show()


@bt.timeshow
def azsnr(opath, npath, target):
    v_light = 299792458
    dtr = pi / 180
    obs, s_type, label = read_ofile(opath)
    snr = choice(obs, s_type, target)

    satlin, sow = all_ofile_satname(snr)
    time = []
    for i in range(0, len(snr)):
        time.append(snr[i][0])
        try:
            del snr[i][0], snr[i][0]
        except:
            pass
    for i in range(0, len(snr)):
        snr[i] = snr[i][0]

    """
    单历元内SNR数据如下
    [
    [42.9, 46.2, 41.77, 0.0, 0.0], 
    [40.81, 47.05, 42.9, 0.0, 0.0], 
    [45.23, 47.07, 43.44, 0.0, 0.0], 
    [40.69, 45.97, 41.61, 0.0, 0.0], 
    [40.68, 42.68, 37.17, 0.0, 0.0], 
    [41.92, 44.77, 42.31, 0.0, 0.0], 
    [47.94, 49.41, 47.33, 0.0, 0.0], 
    [43.0, 48.2, 43.39, 0.0, 0.0], 
    [46.09, 49.1, 46.13, 0.0, 0.0], 
    [43.82, 43.96, 40.73, 0.0, 0.0], 
    [43.82, 0.0, 44.63, 38.17, 46.71], 
    [49.41, 0.0, 49.37, 43.64, 49.83], 
    [51.48, 0.0, 50.28, 46.3, 50.62], 
    [52.63, 0.0, 50.8, 47.28, 51.06], 
    [43.4, 0.0, 42.98, 38.91, 42.78], 
    [43.37, 0.0, 41.59, 38.81, 43.78], 
    [50.08, 0.0, 49.77, 45.3, 50.44], 
    [46.16, 0.0, 47.25, 42.44, 47.97], 
    [44.76, 0.0, 46.54, 41.34, 47.83], 
    [48.52, 0.0, 46.81, 0.0, 0.0], 
    [46.04, 0.0, 46.12, 0.0, 0.0], 
    [45.11, 0.0, 46.27, 0.0, 0.0]
    ]
    """

    for i in range(0, len(snr)):
        for j in range(0, len(snr[i])):
            # snr[i][j]=snr[i][j][0]
            snr[i][j] = max(snr[i][j])
            # for z in range(1, len(snr[i][j])):
            #     del snr[i][j][-1]
    


    Eph = readnfile(npath)
    pos = bt.all_float_list(bt.separate(str(label['APPROX POSITION XYZ'][0])))
    en = len(snr)
    B, L, H = toBLH(pos[0], pos[1], pos[2])
    Transformationmatrix = array([[-sin(B) * cos(L), -sin(B) * sin(L), cos(B)],
                                  [-sin(L), cos(L), 0],
                                  [cos(B) * cos(L), cos(B) * sin(L), sin(B)]
                                  ])
    Surveystation = []
    mid = []
    for q in range(0, en):
        weeksow = sow[q]
        sats = satlin[q]

        a = recpo_lsa(snr[q], sats, weeksow, Eph)

        # a为卫星位置
        mid = []
        to = (180 / pi)
        for i in range(0, len(a)):
            cha = array([[a[i][0] - pos[0]], [a[i][1] - pos[1]], [a[i][2] - pos[2]]])
            tar = dot(Transformationmatrix, cha)
            tar = tar.tolist()
            tar = [tar[0][0], tar[1][0], tar[2][0]]
            # 地球自转改正
            traceltime = sqrt(tar[0] ** 2 + tar[1] ** 2 + tar[2] ** 2) / v_light
            #print(a[i])
            sateloc = e_r_corr(traceltime, a[i])
            cha = array([[sateloc[0] - pos[0]], [sateloc[1] - pos[1]], [sateloc[2] - pos[2]]])
            tar = dot(Transformationmatrix, cha)
            tar = tar.tolist()
            tar = [tar[0][0], tar[1][0], tar[2][0]]
            r = sqrt(tar[0] ** 2 + tar[1] ** 2 + tar[2] ** 2)
            ele = 90 - arcsin(tar[2] / r) * to
            azi = arctan2(tar[1], tar[0])
            if ele > 90:
                mid.append([0, 0, 0])
            else:
                mid.append([azi, ele, snr[q][i]])
        Surveystation.append(mid)
    return time, satlin, Surveystation


@bt.timeshow
def azintsnr(ae, satlin, type):
    """
    Args:
        ae: 卫星高度角方位角 数据来源azsnr
        satlin: 卫星号 数据来源azsnr
        type: 转换类型 1为笛卡尔 -1为极坐标

    Returns:satellite字典 下接drawazsnr

    """
    satellite = {
        "G": {
            "01": [[], [], []], "02": [[], [], []], "03": [[], [], []], "04": [[], [], []], "05": [[], [], []],
            "06": [[], [], []], "07": [[], [], []], "08": [[], [], []], "09": [[], [], []], "10": [[], [], []],
            "11": [[], [], []], "12": [[], [], []], "13": [[], [], []], "14": [[], [], []], "15": [[], [], []],
            "16": [[], [], []], "17": [[], [], []], "18": [[], [], []], "19": [[], [], []], "20": [[], [], []],
            "21": [[], [], []], "22": [[], [], []], "23": [[], [], []], "24": [[], [], []], "25": [[], [], []],
            "26": [[], [], []], "27": [[], [], []], "28": [[], [], []], "29": [[], [], []], "30": [[], [], []],
            "31": [[], [], []], "32": [[], [], []], "33": [[], [], []], "34": [[], [], []], "35": [[], [], []],
            "36": [[], [], []], "37": [[], [], []], "38": [[], [], []], "39": [[], [], []], "40": [[], [], []],
            "41": [[], [], []], "42": [[], [], []], "43": [[], [], []], "44": [[], [], []], "45": [[], [], []],
            "46": [[], [], []], "47": [[], [], []], "48": [[], [], []], "49": [[], [], []], "50": [[], [], []],
            "51": [[], [], []], "52": [[], [], []], "53": [[], [], []], "54": [[], [], []], "55": [[], [], []],
            "56": [[], [], []], "57": [[], [], []], "58": [[], [], []], "59": [[], [], []], "60": [[], [], []], },
        "E": {
            "01": [[], [], []], "02": [[], [], []], "03": [[], [], []], "04": [[], [], []], "05": [[], [], []],
            "06": [[], [], []], "07": [[], [], []], "08": [[], [], []], "09": [[], [], []], "10": [[], [], []],
            "11": [[], [], []], "12": [[], [], []], "13": [[], [], []], "14": [[], [], []], "15": [[], [], []],
            "16": [[], [], []], "17": [[], [], []], "18": [[], [], []], "19": [[], [], []], "20": [[], [], []],
            "21": [[], [], []], "22": [[], [], []], "23": [[], [], []], "24": [[], [], []], "25": [[], [], []],
            "26": [[], [], []], "27": [[], [], []], "28": [[], [], []], "29": [[], [], []], "30": [[], [], []],
            "31": [[], [], []], "32": [[], [], []], "33": [[], [], []], "34": [[], [], []], "35": [[], [], []],
            "36": [[], [], []], "37": [[], [], []], "38": [[], [], []], "39": [[], [], []], "40": [[], [], []],
            "41": [[], [], []], "42": [[], [], []], "43": [[], [], []], "44": [[], [], []], "45": [[], [], []],
            "46": [[], [], []], "47": [[], [], []], "48": [[], [], []], "49": [[], [], []], "50": [[], [], []],
            "51": [[], [], []], "52": [[], [], []], "53": [[], [], []], "54": [[], [], []], "55": [[], [], []],
            "56": [[], [], []], "57": [[], [], []], "58": [[], [], []], "59": [[], [], []], "60": [[], [], []], },
        "C": {
            "01": [[], [], []], "02": [[], [], []], "03": [[], [], []], "04": [[], [], []], "05": [[], [], []],
            "06": [[], [], []], "07": [[], [], []], "08": [[], [], []], "09": [[], [], []], "10": [[], [], []],
            "11": [[], [], []], "12": [[], [], []], "13": [[], [], []], "14": [[], [], []], "15": [[], [], []],
            "16": [[], [], []], "17": [[], [], []], "18": [[], [], []], "19": [[], [], []], "20": [[], [], []],
            "21": [[], [], []], "22": [[], [], []], "23": [[], [], []], "24": [[], [], []], "25": [[], [], []],
            "26": [[], [], []], "27": [[], [], []], "28": [[], [], []], "29": [[], [], []], "30": [[], [], []],
            "31": [[], [], []], "32": [[], [], []], "33": [[], [], []], "34": [[], [], []], "35": [[], [], []],
            "36": [[], [], []], "37": [[], [], []], "38": [[], [], []], "39": [[], [], []], "40": [[], [], []],
            "41": [[], [], []], "42": [[], [], []], "43": [[], [], []], "44": [[], [], []], "45": [[], [], []],
            "46": [[], [], []], "47": [[], [], []], "48": [[], [], []], "49": [[], [], []], "50": [[], [], []],
            "51": [[], [], []], "52": [[], [], []], "53": [[], [], []], "54": [[], [], []], "55": [[], [], []],
            "56": [[], [], []], "57": [[], [], []], "58": [[], [], []], "59": [[], [], []], "60": [[], [], []],
        },
        "R": {
            "01": [[], [], []], "02": [[], [], []], "03": [[], [], []], "04": [[], [], []], "05": [[], [], []],
            "06": [[], [], []], "07": [[], [], []], "08": [[], [], []], "09": [[], [], []], "10": [[], [], []],
            "11": [[], [], []], "12": [[], [], []], "13": [[], [], []], "14": [[], [], []], "15": [[], [], []],
            "16": [[], [], []], "17": [[], [], []], "18": [[], [], []], "19": [[], [], []], "20": [[], [], []],
            "21": [[], [], []], "22": [[], [], []], "23": [[], [], []], "24": [[], [], []], "25": [[], [], []],
            "26": [[], [], []], "27": [[], [], []], "28": [[], [], []], "29": [[], [], []], "30": [[], [], []],
            "31": [[], [], []], "32": [[], [], []], "33": [[], [], []], "34": [[], [], []], "35": [[], [], []],
            "36": [[], [], []], "37": [[], [], []], "38": [[], [], []], "39": [[], [], []], "40": [[], [], []],
            "41": [[], [], []], "42": [[], [], []], "43": [[], [], []], "44": [[], [], []], "45": [[], [], []],
            "46": [[], [], []], "47": [[], [], []], "48": [[], [], []], "49": [[], [], []], "50": [[], [], []],
            "51": [[], [], []], "52": [[], [], []], "53": [[], [], []], "54": [[], [], []], "55": [[], [], []],
            "56": [[], [], []], "57": [[], [], []], "58": [[], [], []], "59": [[], [], []], "60": [[], [], []],
        },
        "J": {
            "01": [[], [], []], "02": [[], [], []], "03": [[], [], []], "04": [[], [], []], "05": [[], [], []],
            "06": [[], [], []], "07": [[], [], []], "08": [[], [], []], "09": [[], [], []], "10": [[], [], []],
            "11": [[], [], []], "12": [[], [], []], "13": [[], [], []], "14": [[], [], []], "15": [[], [], []],
            "16": [[], [], []], "17": [[], [], []], "18": [[], [], []], "19": [[], [], []], "20": [[], [], []],
            "21": [[], [], []], "22": [[], [], []], "23": [[], [], []], "24": [[], [], []], "25": [[], [], []],
            "26": [[], [], []], "27": [[], [], []], "28": [[], [], []], "29": [[], [], []], "30": [[], [], []],
            "31": [[], [], []], "32": [[], [], []], "33": [[], [], []], "34": [[], [], []], "35": [[], [], []],
            "36": [[], [], []], "37": [[], [], []], "38": [[], [], []], "39": [[], [], []], "40": [[], [], []],
            "41": [[], [], []], "42": [[], [], []], "43": [[], [], []], "44": [[], [], []], "45": [[], [], []],
            "46": [[], [], []], "47": [[], [], []], "48": [[], [], []], "49": [[], [], []], "50": [[], [], []],
            "51": [[], [], []], "52": [[], [], []], "53": [[], [], []], "54": [[], [], []], "55": [[], [], []],
            "56": [[], [], []], "57": [[], [], []], "58": [[], [], []], "59": [[], [], []], "60": [[], [], []],
        },
    }
    if type == 1:  # 笛卡尔
        for i in range(0, len(ae)):
            for j in range(0, len(ae[i])):
                m = ae[i][j][0]
                if m < 0:
                    m = m + 2 * math.pi
                satellite[satlin[i][j][0]][satlin[i][j][1:]][0].append(ae[i][j][2])
                satellite[satlin[i][j][0]][satlin[i][j][1:]][1].append(m)
                satellite[satlin[i][j][0]][satlin[i][j][1:]][2].append(90 - ae[i][j][1])
    elif type == -1:  # 极坐标
        for i in range(0, len(ae)):
            for j in range(0, len(ae[i])):
                try:
                    satellite[satlin[i][j][0]][satlin[i][j][1:]][0].append(ae[i][j][2])
                    satellite[satlin[i][j][0]][satlin[i][j][1:]][1].append(ae[i][j][0])
                    satellite[satlin[i][j][0]][satlin[i][j][1:]][2].append(ae[i][j][1])
                except:
                    pass
    return satellite


@bt.timeshow
def drawazsnr(satellite, type):
    """
    Args:
        satellite: 卫星高度角方位角字典，由azintsnr()函数给出
        type: 绘制类型 1为笛卡尔 -1为极坐标

    Returns:绘制卫星天空轨迹图，点的颜色代表信噪比

    """
    if type == 1:
        azi = []
        ele = []
        snr = []
        name = []
        namelocal = []
        for key in satellite.keys():
            for keys in satellite[key].keys():
                if len(satellite[key][keys][0]) > 1 and (key != "C" and keys not in ["56", "57", "58", "59"]):
                    for i in range(0, len(satellite[key][keys][0])):
                        azi.append(satellite[key][keys][1][i])
                        ele.append(satellite[key][keys][2][i])
                        snr.append(satellite[key][keys][0][i])
                    name.append((key, keys))
                    namelocal.append((satellite[key][keys][1][0], satellite[key][keys][2][0]))

        for i in range(0, len(name)):
             plt.annotate(name[i][0] + name[i][1], (namelocal[i][0], namelocal[i][1]))

        c = plt.scatter(azi, ele, c=snr, s=1, cmap='jet')
        c = plt.scatter(azi, ele, c=snr, s=1, cmap='jet')
        # c = plt.hexbin(azi, ele,c=snr, cmap='jet')
        cb = plt.colorbar(c, )
        plt.xlabel("azimuth")
        plt.ylabel("elevation")
        plt.xlim(0, math.pi * 2)
        plt.ylim(0, 90)
        xticks = [0, 1 * math.pi / 3, 2 * math.pi / 3, math.pi, 4 * math.pi / 3, 5 * math.pi / 3, 2 * math.pi]
        xticklabel = ["0", "60", "120", "180", "240", "300", "360"]
        plt.xticks(xticks, xticklabel)
        plt.title("SNR value")
        plt.legend()
        plt.savefig("GNSS sky map.png", dpi=1000)
        plt.show()
    elif type == -1:
        ax = plt.subplot(111, projection='polar')
        ax.set_theta_direction(-1)
        ax.set_theta_zero_location('N')
        ax.yaxis.set_label_position('right')
        ax.tick_params('y', labelleft=False)
        ax.grid(linestyle='--')
        labels = ['N', '45°', 'E', '135°', 'S', '225°', 'W', '315°']
        ax.set_thetagrids(range(0, 360, 45), labels, fontweight='semibold')
        azi, ele, snr, name, namelocal = [], [], [], [], []
        bds3 = ["19","20","23","24","25","26","27","29","30","32","34","35","36","38","39",
                "40","43","44","45","46"]

        for key in satellite.keys():
            for keys in satellite[key].keys():
                if keys in bds3:
                    if len(satellite[key][keys][0]) > 1:
                        azi = azi + satellite[key][keys][1]
                        ele = ele + satellite[key][keys][2]
                        snr = snr + satellite[key][keys][0]
                        name.append((key, keys))
                        namelocal.append((satellite[key][keys][1][0], satellite[key][keys][2][0]))
                # if len(satellite[key][keys][0]) > 1:
                #     azi = azi + satellite[key][keys][1]
                #     ele = ele + satellite[key][keys][2]
                #     snr = snr + satellite[key][keys][0]
                #     name.append((key, keys))
                #     namelocal.append((satellite[key][keys][1][0], satellite[key][keys][2][0]))
        for i in range(0, len(name)):
            plt.annotate(name[i][0] + name[i][1], (namelocal[i][0], namelocal[i][1]),)
        # demoazi,demoele,demosnr = [],[],[]
        # for i in range(0,len(snr)):
        #     if snr[i]<=30:
        #         demosnr.append(snr[i])
        #         demoazi.append(azi[i])
        #         demoele.append(ele[i])

        #c= plt.scatter(demoazi,demoele,c = demosnr,s=1,cmap='rainbow')

        azi.append(-1)
        ele.append(90)
        snr.append(0)


        c= plt.scatter(azi,ele,c = snr,s=1,cmap='rainbow',label="Local_SNR")
        cb = plt.colorbar(c, )
        SATAZ = [0, 0, 0, 0, 0, 0]
        SATEL = [0, 15, 30, 45, 60, 75]
        SATID = ['90', '75', '60', '45', '30', '15']
        for i in range(0, 6):
            ax.text(SATAZ[i], SATEL[i], SATID[i])
        ax.set_rticks(range(0, 105, 15))
        #plt.title("HF-2023-remote-4-G")
        plt.legend()
        plt.show()


@bt.timeshow
def density_calc(x, y):
    res = empty((360, 90), dtype=float32)
    for i in range(0, len(x)):
        try:
            res[math.floor(x[i])][math.floor(y[i])] = res[math.floor(x[i])][math.floor(y[i])] + 1
        except:
            pass
    del x, y
    num = 0
    res2 = empty((360, 90))
    for i in range(5, 355):
        for j in range(5, 85):
            for m in [-5, -4, -3, -2, -1, 0, 1, 1, 2, 3, 4, 5]:
                for n in [-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5]:
                    if res[i + m][j + n] >= 1:
                        num = num + 1
            res2[i][j] = res2[i][j] + num
            num = 0
    hot = []
    azi = []
    ele = []
    for i in range(0, len(res2)):
        for j in range(0, len(res2[i])):
            if abs(res2[i][j]) < 1:
                hot.append(res2[i][j])
                azi.append(i)
                ele.append(j)
    del res
    return hot, azi, ele


@bt.timeshow
def drawnumhot(satellite, type):
    if type == 1:
        azi = []
        ele = []
        name = []
        namelocal = []
        rou = 180 / pi
        for key in satellite.keys():
            for keys in satellite[key].keys():
                if len(satellite[key][keys][0]) > 1 and (key != "C" and keys not in ["56", "57", "58", "59"]):
                    for i in range(0, len(satellite[key][keys][0])):
                        azi.append(satellite[key][keys][1][i] * rou)
                    ele = ele + satellite[key][keys][2]
                    name.append((key, keys))
                    namelocal.append((satellite[key][keys][1][0], satellite[key][keys][2][0]))
        for i in range(0, len(name)):
            plt.annotate(name[i][0] + name[i][1], (namelocal[i][0] * rou, namelocal[i][1]))

        r, a, e = density_calc(azi, ele)
        colormap = plt.get_cmap("jet")
        # plt.scatter(a, e, c=r, cmap=colormap, norm=colors.LogNorm(vmin=min(r), vmax=max(r)))
        plt.scatter(a, e, c=r, cmap=colormap, )
        plt.grid(linestyle='--', color="grey")
        plt.xlabel("azimuth")
        plt.ylabel("elevation")
        plt.xlim(0, 360)
        plt.ylim(0, 90)
        plt.legend()
        plt.savefig("GNSS sky map.png", dpi=1000)
        plt.show()
    elif type == -1:
        ax = plt.subplot(111, projection='polar')
        ax.set_theta_direction(-1)
        ax.set_theta_zero_location('N')
        ax.yaxis.set_label_position('right')
        ax.tick_params('y', labelleft=False)
        ax.grid(linestyle='--')
        labels = ['N', '45°', 'E', '135°', 'S', '225°', 'W', '315°']
        ax.set_thetagrids(range(0, 360, 45), labels, fontweight='semibold')
        azi, ele, snr, name, namelocal = [], [], [], [], []
        for key in satellite.keys():
            for keys in satellite[key].keys():
                if len(satellite[key][keys][0]) > 1:
                    azi = azi + satellite[key][keys][1]
                    ele = ele + satellite[key][keys][2]
                    snr = snr + satellite[key][keys][0]
                    name.append((key, keys))
                    namelocal.append((satellite[key][keys][1][0], satellite[key][keys][2][0]))
        for i in range(0, len(name)):
            plt.annotate(name[i][0] + name[i][1], (namelocal[i][0], namelocal[i][1]))
        # c= plt.scatter(azi,ele,c = snr,s=1,cmap='jet')
        c = plt.hexbin(azi, ele, )
        cb = plt.colorbar(c, )
        SATAZ = [0, 0, 0, 0, 0, 0]
        SATEL = [0, 15, 30, 45, 60, 75]
        SATID = ['90', '75', '60', '45', '30', '15']
        for i in range(0, 6):
            ax.text(SATAZ[i], SATEL[i], SATID[i])
        ax.set_rticks(range(0, 105, 15))
        plt.legend()
        plt.show()


@bt.timeshow
def azlocal(opath, npath, target):
    v_light = 299792458
    dtr = pi / 180
    obs, s_type, label = read_ofile(opath)
    obs = choice(obs, s_type, target)
    satlin, sow = all_ofile_satname(obs)
    time = []
    for i in range(0, len(obs)):
        time.append(obs[i][0])
        try:
            del obs[i][0], obs[i][0]
        except:
            pass
    for i in range(0, len(obs)):
        obs[i] = obs[i][0]
    for i in range(0, len(obs)):
        for j in range(0, len(obs[i])):
            for z in range(1, len(obs[i][j])):
                del obs[i][j][-1]
    Eph = readnfile(npath)
    pos = bt.all_float_list(bt.separate(str(label['APPROX POSITION XYZ'][0])))
    en = len(obs)
    # B, L, H = toBLH(pos[0], pos[1], pos[2])
    B, L, H = 90, 0, 0
    Transformationmatrix = array([[-sin(B) * cos(L), -sin(B) * sin(L), cos(B)],
                                  [-sin(L), cos(L), 0],
                                  [cos(B) * cos(L), cos(B) * sin(L), sin(B)]
                                  ])
    Surveystation = []
    mid = []
    for q in range(0, en):
        weeksow = sow[q]
        sats = satlin[q]
        a = recpo_lsa(obs[q], sats, weeksow, Eph)
        # a为卫星位置
        mid = []
        to = (180 / pi)
        for i in range(0, len(a)):
            cha = array([[a[i][0] - pos[0]], [a[i][1] - pos[1]], [a[i][2] - pos[2]]])
            tar = dot(Transformationmatrix, cha)
            tar = tar.tolist()
            tar = [tar[0][0], tar[1][0], tar[2][0]]
            # 地球自转改正
            traceltime = sqrt(tar[0] ** 2 + tar[1] ** 2 + tar[2] ** 2) / v_light
            sateloc = e_r_corr(traceltime, a[i])
            cha = array([[sateloc[0] - pos[0]], [sateloc[1] - pos[1]], [sateloc[2] - pos[2]]])
            tar = dot(Transformationmatrix, cha)
            tar = tar.tolist()
            tar = [tar[0][0], tar[1][0], tar[2][0]]
            r = sqrt(tar[0] ** 2 + tar[1] ** 2 + tar[2] ** 2)
            ele = 90 - arcsin(tar[2] / r) * to
            azi = arctan2(tar[1], tar[0])
            if ele > 90:
                mid.append([0, 0])
            else:
                mid.append([azi, ele, ])
        Surveystation.append(mid)
    return time, satlin, Surveystation


def fake_recpo(sats, time, Eph):
    col_Eph = []
    m = len(sats)
    for t in range(0, m):
        col = find_eph(Eph, sats[t], time)
        col_Eph.append(col)
    xpos = []
    n = len(col_Eph)
    if n == 0:
        return nan
    for i in range(0, n):
        k = col_Eph[i]
        if sats[i][0] == "R":
            a = check_t(time)
            if a < 0:
                a = 86400 + a
            h = check_t(a - Eph[k]['TEO']) + 3 * 3600 - 15 * 60
            X = gloposition(Eph[k]['x'], Eph[k]['y'], Eph[k]['z'], Eph[k]['xv'], Eph[k]['yv'], Eph[k]['zv'],
                            Eph[k]['xa'], Eph[k]['ya'], Eph[k]['za'], h)
            xpos.append(X)
        else:
            tx_RAW = time
            toc = Eph[k]['TEO']
            dt = check_t(tx_RAW - toc)
            try:
                tcorr = (Eph[k]['卫星钟漂速度参数（s/s/s）'] * dt + Eph[k]['卫星钟漂（s/s）']) * dt + Eph[k]['卫星钟差（s）']
            except:
                tcorr = 0
            tx_GPS = tx_RAW - tcorr
            dt = check_t(tx_GPS - toc)
            try:
                tcorr = (Eph[k]['卫星钟漂速度参数（s/s/s）'] * dt + Eph[k]['卫星钟漂（s/s）']) * dt + Eph[k]['卫星钟差（s）']
            except:
                tcorr = 0
            tx_GPS = tx_RAW - tcorr
            if sats[i] in ["C01", "C02", "C03", "C04", "C05", "C59", "C60"]:
                X = satposgeo(tx_GPS, Eph[k])
            else:
                X = satpos(tx_GPS, Eph[k])
            xpos.append(X)
    return xpos


@bt.timeshow
def fake_skymap(nfile):
    v_light = 299792458
    eph = readnfile(nfile)
    geph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "G"]
    ceph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "C"]
    eeph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "E"]
    reph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "R"]
    jeph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "J"]
    ieph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "I"]
    seph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "S"]

    co_g = {"G": {"GM": 3.986004418e14, "Omegae_dot": 7.2921151467e-5},
            "C": {"GM": 3.986004418e14, "Omegae_dot": 7.2921150000e-5},
            "R": {"GM": 3.986004400e14, "Omegae_dot": 7.2921150000e-5},
            "E": {"GM": 3.986004418e14, "Omegae_dot": 7.2921151467e-5},
            "J": {"GM": 3.986004418e14, "Omegae_dot": 7.2921151467e-5},
            }

    s_eph_name = list(set([eph[i]['卫星PRN号'] for i in range(0, len(eph))]))
    gsate = sort([i for i in s_eph_name if i[0] == "G"])
    csate = sort([i for i in s_eph_name if i[0] == "C"])
    esate = sort([i for i in s_eph_name if i[0] == "E"])
    rsate = sort([i for i in s_eph_name if i[0] == "R"])
    jsate = sort([i for i in s_eph_name if i[0] == "J"])
    isate = sort([i for i in s_eph_name if i[0] == "I"])
    ssate = sort([i for i in s_eph_name if i[0] == "S"])
    del s_eph_name, eph
    gtimestart = gpstime(julday(bt.all_float_list(bt.separate(geph[0]['历元']))))
    gtimeend = gpstime(julday(bt.all_float_list(bt.separate(geph[-1]['历元']))))
    timelength = gtimeend - gtimestart

    gobs = [[float(gtimestart) + float(i)] for i in range(0, int(timelength))]
    en, mid, dtr = len(gobs), [], 180 / pi
    B, L, H = 0, 0, 0
    Transformationmatrix = array([[-sin(B) * cos(L), -sin(B) * sin(L), cos(B)],
                                  [-sin(L), cos(L), 0],
                                  [cos(B) * cos(L), cos(B) * sin(L), sin(B)]
                                  ])
    #pos = Blhtoxyz84(B,L,H)
    pos = [6378137,0,0]
    Surveystation = []
    satlin = [gsate for i in range(len(gobs))]
    time = [gobs[q][0] for q in range(0, len(gobs))]

    for q in range(0, en):
        weeksow, sats = time[q], satlin[q]
        posline = fake_recpo(sats, weeksow, geph)
        mid =[]
        for i in range(0, len(posline)):
            tar = array([[posline[i][0] - pos[0]], [posline[i][1] - pos[1]], [posline[i][2] - pos[2]]])
            # tar = dot(Transformationmatrix, cha).tolist()
            # tar = [tar[0][0], tar[1][0], tar[2][0]]
            # 地球自转改正
            # traceltime = sqrt(tar[0] ** 2 + tar[1] ** 2 + tar[2] ** 2) / v_light
            # sateloc = e_r_corr(traceltime, posline[i])
            # cha = array([[sateloc[0] - pos[0]], [sateloc[1] - pos[1]], [sateloc[2] - pos[2]]])
            # tar = dot(Transformationmatrix, cha).tolist()
            # tar = [tar[0][0], tar[1][0], tar[2][0]]
            # tar =array([[posline[i][0] - pos[0]], [posline[i][1] - pos[1]], [posline[i][2] - pos[2]]])
            r = sqrt(tar[0] ** 2 + tar[1] ** 2 + tar[2] ** 2)
            ele = 90 - arcsin(tar[2] / r) * dtr
            azi = arctan2(tar[1], tar[0])
            if ele > 90:
                mid.append([0, 0])
            else:
                mid.append([azi, ele])
        Surveystation.append(mid)
    return time, satlin, Surveystation


@bt.timeshow
def map3d(nfile):
    """
    Args:
        nfile: 混合星历文件

    Returns:绘制GPS卫星轨道三维图
    """
    v_light = 299792458
    eph = readnfile(nfile)
    geph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "G"]
    ceph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "C"]
    eeph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "E"]
    # reph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "R"]
    jeph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "J"]
    ieph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "I"]
    seph = [eph[i] for i in range(0, len(eph)) if eph[i]['卫星PRN号'][0] == "S"]

    co_g = {"G": {"GM": 3.986004418e14, "Omegae_dot": 7.2921151467e-5},
            "C": {"GM": 3.986004418e14, "Omegae_dot": 7.2921150000e-5},
            # "R": {"GM": 3.986004400e14, "Omegae_dot": 7.2921150000e-5},
            "E": {"GM": 3.986004418e14, "Omegae_dot": 7.2921151467e-5},
            "J": {"GM": 3.986004418e14, "Omegae_dot": 7.2921151467e-5},
            }

    s_eph_name = list(set([eph[i]['卫星PRN号'] for i in range(0, len(eph))]))
    gsate = sort([i for i in s_eph_name if i[0] == "G"])
    csate = sort([i for i in s_eph_name if i[0] == "C"])
    esate = sort([i for i in s_eph_name if i[0] == "E"])
    # rsate = sort([i for i in s_eph_name if i[0] == "R"])
    jsate = sort([i for i in s_eph_name if i[0] == "J"])
    isate = sort([i for i in s_eph_name if i[0] == "I"])
    ssate = sort([i for i in s_eph_name if i[0] == "S"])
    del s_eph_name, eph
    gtimestart = gpstime(julday(bt.all_float_list(bt.separate(geph[0]['历元']))))
    gtimeend = gpstime(julday(bt.all_float_list(bt.separate(geph[-1]['历元']))))
    timelength = gtimeend - gtimestart

    gobs = [[float(gtimestart) + float(i)] for i in range(0, int(timelength))]
    en, mid, dtr = len(gobs), [], 180 / pi
    B, L, H = 60, 0, 0
    Transformationmatrix = array([[-sin(B) * cos(L), -sin(B) * sin(L), cos(B)],
                                  [-sin(L), cos(L), 0],
                                  [cos(B) * cos(L), cos(B) * sin(L), sin(B)]
                                  ])
    # pos = Blhtoxyz84(B,L,H)
    pos = [-2339499.7543,4555687.5821,3789167.9972]
    Surveystation = []
    satlin = [gsate for i in range(len(gobs))]
    time = [gobs[q][0] for q in range(0, len(gobs))]
    x,y,z=[],[],[]
    for q in range(0, en):
        weeksow, sats = time[q], satlin[q]
        posline = fake_recpo(sats, weeksow, geph)
        if q%10 <0.5:
            for i in range(len(posline)):
                x.append(posline[i][0])
                y.append(posline[i][1])
                z.append(posline[i][2])
    ax =plt.subplot(projection='3d')
    ax.scatter(x,y,z,s=0.1)
    plt.show()

def convert_time(time_list:list):
    """
    将历元时间转换成Datatime
    """
    from datetime import datetime

    datetime_obj = datetime(year=int(time_list[0]), month=int(time_list[1]), day=int(time_list[2]),
                        hour=int(time_list[3]), minute=int(time_list[4]), second=int(time_list[5]),
                        microsecond=int((time_list[5]*10)%10)*100000)


    return datetime_obj

def Data_completion_rate(file_path:str):
    """
    数据完整率
    
    """
    import numpy as np
    cal,s_type,labels = read_ofile(file_path)
    satellite_set = list(set([i[j][0] for i in cal for j in range(1, len(i))]))

    satellites_dict = {}
    for i in labels['SYS / # / OBS TYPES']:
        satellites_dict.update({i[0]:{}})
    
    
    GNSS_obs_type = {}
    for i in labels['SYS / # / OBS TYPES']:
        sp = i.split(" ")
        GNSS_obs_type.update({sp[0]: sp[5:]})
    for i in satellite_set:
        obs_type = i[0]
        satellites_dict[i[0]].update({i: {}})
        for j in range(0, len(GNSS_obs_type[obs_type])):
            satellites_dict[i[0]][i].update({GNSS_obs_type[obs_type][j]: {"time": [], "value": []}})
    for i in cal:
        time = convert_time(i[0])
        for j in range(1, len(i)):
            mid_sate = i[j][0]
            mid_type = mid_sate[0]
            content = i[j][1:]
            for k in range(0, len(content)):
                satellites_dict[mid_type][mid_sate][GNSS_obs_type[mid_type][k]]["time"].append(time)
                satellites_dict[mid_type][mid_sate][GNSS_obs_type[mid_type][k]]["value"].append(content[k])
    syst = {}
    res_put={}

    for i in GNSS_obs_type.keys():
        syst.update({i: {}})
        res_put.update({i:{}})
        for j in GNSS_obs_type[i]:
            syst[i].update({j: []})
            res_put[i].update({j:0})
            
    
    for i in satellites_dict.keys():

        for j in satellites_dict[i].keys():
            for k in satellites_dict[i][j].keys():
                syst[i][k].extend(satellites_dict[i][j][k]['value'])

    for i in syst.keys():
        for j in syst[i].keys():
            a = syst[i][j]
            b = sum(np.where(np.array(a) > 1, 1, 0)) / len(a)
            # print(i, j, "的数据完整率为：", b * 100, "%")
            res_put[i][j]=round(b*100,3)
    return res_put,syst
# ofile = r"F:\静态文件\2022.3.4逸夫\2022.3.5逸夫楼\Rinex\3392318064B4.22o"
# nfile = r"F:\静态文件\2022.3.4逸夫\2022.3.5逸夫楼\Rinex\3392318064B4.22p"
#
# time, satlin, Surveystation=azsnr(ofile,nfile,(["G","C","E","J"],["S"]))
# s=azintsnr(Surveystation,satlin,-1)
# drawazsnr(s,-1)

